home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / fsolve.cc < prev    next >
C/C++ Source or Header  |  1997-05-26  |  7KB  |  329 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include <string>
  28.  
  29. #include <iostream.h>
  30.  
  31. #include "NLEqn.h"
  32.  
  33. #include "defun-dld.h"
  34. #include "error.h"
  35. #include "gripes.h"
  36. #include "help.h"
  37. #include "pager.h"
  38. #include "pt-fvc.h"
  39. #include "oct-obj.h"
  40. #include "utils.h"
  41. #include "variables.h"
  42.  
  43. // Global pointer for user defined function required by hybrd1.
  44. static tree_fvc *fsolve_fcn;
  45.  
  46. static NLEqn_options fsolve_opts;
  47.  
  48. int
  49. hybrd_info_to_fsolve_info (int info)
  50. {
  51.   switch (info)
  52.     {
  53.     case -1:
  54.       info = -2;
  55.       break;
  56.  
  57.     case 0:
  58.       info = -1;
  59.       break;
  60.  
  61.     case 1:
  62.       break;
  63.  
  64.     case 2:
  65.       info = 4;
  66.       break;
  67.  
  68.     case 3:
  69.     case 4:
  70.     case 5:
  71.       info = 3;
  72.       break;
  73.  
  74.     default:
  75.       panic_impossible ();
  76.       break;
  77.     }
  78.  
  79.   return info;
  80. }
  81.  
  82. ColumnVector
  83. fsolve_user_function (const ColumnVector& x)
  84. {
  85.   ColumnVector retval;
  86.  
  87.   int n = x.capacity ();
  88.  
  89.   octave_value_list args;
  90.   args.resize (1);
  91.  
  92.   if (n > 1)
  93.     {
  94.       Matrix m (n, 1);
  95.       for (int i = 0; i < n; i++)
  96.     m (i, 0) = x (i);
  97.       octave_value vars (m);
  98.       args(0) = vars;
  99.     }
  100.   else
  101.     {
  102.       double d = x (0);
  103.       octave_value vars (d);
  104.       args(0) = vars;
  105.     }
  106.  
  107.   if (fsolve_fcn)
  108.     {
  109.       octave_value_list tmp = fsolve_fcn->eval (0, 1, args);
  110.       if (tmp.length () > 0 && tmp(0).is_defined ())
  111.     {
  112.       retval = tmp(0).vector_value ();
  113.  
  114.       if (error_state || retval.length () <= 0)
  115.         gripe_user_supplied_eval ("fsolve");
  116.     }
  117.       else
  118.     gripe_user_supplied_eval ("fsolve");
  119.     }
  120.  
  121.   return retval;
  122. }
  123.  
  124. DEFUN_DLD (fsolve, args, nargout,
  125.   "Solve nonlinear equations using Minpack.  Usage:\n\
  126. \n\
  127.   [X, INFO] = fsolve (F, X0)\n\
  128. \n\
  129. Where the first argument is the name of the  function to call to\n\
  130. compute the vector of function values.  It must have the form\n\
  131. \n\
  132.   y = f (x)\n\
  133. \n\
  134. where y and x are vectors.")
  135. {
  136.   octave_value_list retval;
  137.  
  138.   int nargin = args.length ();
  139.  
  140.   if (nargin != 2 || nargout > 3)
  141.     {
  142.       print_usage ("fsolve");
  143.       return retval;
  144.     }
  145.  
  146.   fsolve_fcn = is_valid_function (args(0), "fsolve", 1);
  147.   if (! fsolve_fcn)
  148.     return retval;
  149.  
  150.   ColumnVector x = args(1).vector_value ();
  151.  
  152.   if (error_state)
  153.     {
  154.       error ("fsolve: expecting vector as second argument");
  155.       return retval;
  156.     }
  157.  
  158.   if (nargin > 2)
  159.     warning ("fsolve: ignoring extra arguments");
  160.  
  161.   if (nargout > 2)
  162.     warning ("fsolve: can't compute path output yet");
  163.  
  164.   NLFunc foo_fcn (fsolve_user_function);
  165.   NLEqn foo (x, foo_fcn);
  166.   foo.set_options (fsolve_opts);
  167.  
  168.   int info;
  169.   ColumnVector soln = foo.solve (info);
  170.  
  171.   info = hybrd_info_to_fsolve_info (info);
  172.  
  173.   retval.resize (nargout ? nargout : 1);
  174.   retval(0) = soln, 1;
  175.  
  176.   if (nargout > 1)
  177.     retval(1) = (double) info;
  178.  
  179.   return retval;
  180. }
  181.  
  182. typedef void (NLEqn_options::*d_set_opt_mf) (double);
  183. typedef double (NLEqn_options::*d_get_opt_mf) (void);
  184.  
  185. #define MAX_TOKENS 1
  186.  
  187. struct NLEQN_OPTIONS
  188. {
  189.   const char *keyword;
  190.   const char *kw_tok[MAX_TOKENS + 1];
  191.   int min_len[MAX_TOKENS + 1];
  192.   int min_toks_to_match;
  193.   d_set_opt_mf d_set_fcn;
  194.   d_get_opt_mf d_get_fcn;
  195. };
  196.  
  197. static NLEQN_OPTIONS fsolve_option_table [] =
  198. {
  199.   { "tolerance",
  200.     { "tolerance", 0, },
  201.     { 1, 0, }, 1,
  202.     NLEqn_options::set_tolerance,
  203.     NLEqn_options::tolerance, },
  204.  
  205.   { 0,
  206.     { 0, 0, },
  207.     { 0, 0, }, 0,
  208.     0, 0, },
  209. };
  210.  
  211. static void
  212. print_fsolve_option_list (ostream& os)
  213. {
  214.   print_usage ("fsolve_options", 1);
  215.  
  216.   os << "\n"
  217.      << "Options for fsolve include:\n\n"
  218.      << "  keyword                                  value\n"
  219.      << "  -------                                  -----\n\n";
  220.  
  221.   NLEQN_OPTIONS *list = fsolve_option_table;
  222.  
  223.   const char *keyword;
  224.   while ((keyword = list->keyword) != 0)
  225.     {
  226.       os.form ("  %-40s ", keyword);
  227.  
  228.       double val = (fsolve_opts.*list->d_get_fcn) ();
  229.       if (val < 0.0)
  230.     os << "computed automatically";
  231.       else
  232.     os << val;
  233.  
  234.       os << "\n";
  235.       list++;
  236.     }
  237.  
  238.   os << "\n";
  239. }
  240.  
  241. static void
  242. set_fsolve_option (const string& keyword, double val)
  243. {
  244.   NLEQN_OPTIONS *list = fsolve_option_table;
  245.  
  246.   while (list->keyword != 0)
  247.     {
  248.       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
  249.                 list->min_toks_to_match, MAX_TOKENS))
  250.     {
  251.       (fsolve_opts.*list->d_set_fcn) (val);
  252.  
  253.       return;
  254.     }
  255.       list++;
  256.     }
  257.  
  258.   warning ("fsolve_options: no match for `%s'", keyword.c_str ());
  259. }
  260.  
  261. static octave_value_list
  262. show_fsolve_option (const string& keyword)
  263. {
  264.   octave_value_list retval;
  265.  
  266.   NLEQN_OPTIONS *list = fsolve_option_table;
  267.  
  268.   while (list->keyword != 0)
  269.     {
  270.       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
  271.                 list->min_toks_to_match, MAX_TOKENS))
  272.     {
  273.       return (fsolve_opts.*list->d_get_fcn) ();
  274.     }
  275.       list++;
  276.     }
  277.  
  278.   warning ("fsolve_options: no match for `%s'", keyword.c_str ());
  279.  
  280.   return retval;
  281. }
  282.  
  283. DEFUN_DLD (fsolve_options, args, ,
  284.   "fsolve_options (KEYWORD, VALUE)\n\
  285. \n\
  286. Set or show options for fsolve.  Keywords may be abbreviated\n\
  287. to the shortest match.")
  288. {
  289.   octave_value_list retval;
  290.  
  291.   int nargin = args.length ();
  292.  
  293.   if (nargin == 0)
  294.     {
  295.       print_fsolve_option_list (octave_stdout);
  296.       return retval;
  297.     }
  298.   else if (nargin == 1 || nargin == 2)
  299.     {
  300.       string keyword = args(0).string_value ();
  301.  
  302.       if (! error_state)
  303.     {
  304.       if (nargin == 1)
  305.         return show_fsolve_option (keyword);
  306.       else
  307.         {
  308.           double val = args(1).double_value ();
  309.  
  310.           if (! error_state)
  311.         {
  312.           set_fsolve_option (keyword, val);
  313.           return retval;
  314.         }
  315.         }
  316.     }
  317.     }
  318.  
  319.   print_usage ("fsolve_options");
  320.  
  321.   return retval;
  322. }
  323.  
  324. /*
  325. ;;; Local Variables: ***
  326. ;;; mode: C++ ***
  327. ;;; End: ***
  328. */
  329.